subroutine hydrocarbon_init()
use molecule
implicit none
   !-------------------------------------------------------------------------
   ! Hydrocarbon chain molecule generator
   !-------------------------------------------------------------------------
   character (len=4) :: HCHC_e(1), CTCT_e(6), CMCM_e(4), CTHC_e(3)
   real(q)           :: HCHC_p(3, 3), CTCT_p(3, 7), CMCM_p(3, 5), CTHC_p(3, 4)
   real(q)           :: HCHC_q(2), CTCT_q(6), CMCM_q(4), CTHC_q(3)
   integer           :: HCHC_n(4, 1), CTCT_n(4, 6), CMCM_n(4, 4), CTHC_n(4, 3)
   integer           :: HCHC_R(2), CTCT_R(2), CMCM_R(2), CTHC_R(2)
   integer           :: HCHC_L(2), CTCT_L(2), CMCM_L(2), CTHC_L(2)
   integer           :: nx = 0, ny = 0, nz = 0, m = 0, n = 0, nq = 1
   integer           :: id, i, j, ii, jj, href
   integer           :: dL, dR, nL, nR
   integer           :: first_CM, first_CT, last_CM, last_CT, flag, ioerr
   real(q)           :: refpos(3)
   character(len=255):: oneline
   !-------------------------------------------------------------------------
   HCHC_e      = (/ "HC" /)
   CTCT_e      = (/ "CT", "HC", "HC", "CT", "HC", "HC" /)
   CMCM_e      = (/ "CM", "HC", "CM", "HC" /)
   CTHC_e      = (/ "CT", "HC", "HC" /)
   HCHC_q      = (/ 0.06D0, 0.115D0/)
   CTCT_q      = (/ -0.12D0, 0.06D0, 0.06D0, -0.12D0, 0.06D0, 0.06D0 /)
   CMCM_q      = (/ -0.115D0, 0.115D0, -0.115D0, 0.115D0 /)
   CTHC_q      = (/ -0.12D0, 0.06D0, 0.06D0 /)

   HCHC_n      = reshape( (/ -1, 0, 0, 0 /), (/ 4, 1 /) )
   CTCT_n      = reshape( (/ 2, 3, 4, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 5, 6, -1, 4, 0, 0, 0, 4, 0, 0, 0 /), (/4, 6/) )
   CMCM_n      = reshape( (/ 2, 3, -1, 0, 1, 0, 0, 0, 1, 4, -1, 0, 3, 0, 0, 0 /), (/4, 4/) )
   CTHC_n      = reshape( (/ 2, 3, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0 /), (/4, 3/) )

   HCHC_R      = (/ 1, 1 /)
   CTCT_R      = (/ 4, 4 /)
   CMCM_R      = (/ 3, 3 /)
   CTHC_R      = (/ 1, 4 /)

   HCHC_L      = (/ 1, 1 /)
   CTCT_L      = (/ 1, 4 /)
   CMCM_L      = (/ 1, 3 /)
   CTHC_L      = (/ 1, 3 /)

   HCHC_p      = reshape( (/ -0.5694D0, 0.D0, 0.9312D0, 0.D0, 0.D0, 0.D0, -0.3D0, 0.D0, 0.4D0 /), (/3, 3/) )
   CTCT_p      = reshape( (/ 0.D0, 0.D0, 0.D0, -0.056D0, -0.8825D0, -0.6694D0, -0.056D0, 0.8825D0, -0.6694D0,       &
                             1.2815D0, 0.D0, 0.8570D0, 1.2816D0, -0.8825D0, 1.5270D0, 1.2916D0, 0.8825D0, 1.5270D0, &
                             2.5603D0, 0.D0, 0.D0 /), (/3, 7/) )
   CMCM_p      = reshape( (/ 0.D0, 0.D0, 0.D0, -0.5694D0, 0.D0, -0.9312D0, 1.334D0, 0.D0, 0.D0, 1.8524D0, 0.D0, 0.9621D0, &
                             2.1487D0, 0.D0, -1.2756D0 /), (/3, 5/))
   CTHC_p      = reshape( (/ 0.D0, 0.D0, 0.D0, -0.056D0, -0.8915D0, -0.6530D0, -0.056D0, 0.8915D0, -0.6694D0,       &
                             0.9117D0, 0.D0, 0.6260D0 /), (/3, 4/) )
   !-------------------------------------------------------------------------
   write(*, '(10X, "We are going to create a molecule of " )')
   write(*, '(15X, " H-(CH2CH2)x-[(CH=CH)m-(CH2CH2)n-]y-(CH=CH)z-(CH2)q-H")')
   write(*, '(10X, "where x, y, z, m, n could be non-negative integers, q in [0, 1].")')
   write(*, '(10X, "Please input the desired x, m, n, y, z, and q [0,0,0,0,0,1]: ", $)')
   read(*, '(A)', iostat=ioerr) oneline
   if (ioerr.eq.0.and.oneline.ne.'') read(oneline, *, iostat=ioerr) nx, m, n, ny, nz, nq
   if (nx.lt.0) nx = 0
   if ( m.lt.0) m  = 0
   if ( n.lt.0) n  = 0
   if (ny.lt.0) ny = 0
   if (nz.lt.0) nz = 0
   if (nq.ne.0) nq = 1
   write(*, '(10X, "The accepted numbers are  :", 6(2X, I0) )') nx, m, n, ny, nz, nq
  
   natom    = 2 + 6*nx + (4*m + 6*n)*ny + 4 * nz + 3*nq
   oneline = "H"
   if (nx == 1) write(oneline, '(A,"-CH2CH2" )') trim(oneline)
   if (nx >  1) write(oneline, '(A,"-(CH2CH2)", I0)') trim(oneline), nx
   if (ny > 0) then
      if (m * n > 0) then
         if (ny > 1) write(oneline, '(A, "-[")' ) trim(oneline)
         if (m == 1) write(oneline, '(A,"-CH=CH")') trim(oneline)
         if (m >  1) write(oneline, '(A,"-(CH=CH)", I0)') trim(oneline),  m
         if (ny == 1) then
            if (n == 1) write(oneline, '(A,"-CH2CH2" )') trim(oneline)
            if (n >  1) write(oneline, '(A,"-(CH2CH2)", I0)') trim(oneline), n
         else
            if (n == 1) write(oneline, '(A,"-CH2CH2-]", I0)') trim(oneline), ny
            if (n >  1) write(oneline, '(A,"-(CH2CH2)", I0, "-]", I0)') trim(oneline), n, ny
         endif
      else
         if (m*ny == 1) write(oneline, '(A,"-CH=CH")') trim(oneline)
         if (m*ny >  1) write(oneline, '(A,"-(CH=CH)", I0)') trim(oneline), m*ny
         if (n*ny == 1) write(oneline, '(A,"-CH2CH2", I0)') trim(oneline)
         if (n*ny >  0) write(oneline, '(A,"-(CH2CH2)", I0)') trim(oneline), n*ny
      endif
   endif
   if (nz == 1) write(oneline, '(A,"-CH=CH")') trim(oneline)
   if (nz >  1) write(oneline, '(A,"-(CH=CH)", I0)') trim(oneline), nz
   if (nq > 0) then
      write(oneline, '(A,"-CH3")') trim(oneline)
   else
      write(oneline, '(A,"-H")') trim(oneline)
   endif
   write(*, '(10X, "The molecule formula reads: ", A, " with ", I0, " atoms.")') trim(oneline), natom
   !
   neighmax = 4
   molname = "hydrocarbon"
   if ( allocated(neigh   )) deallocate( neigh )
   if ( allocated(element )) deallocate( element )
   if ( allocated(molpos  )) deallocate( molpos  )
   if ( allocated(chg     )) deallocate( chg     )
   allocate( element(natom), neigh(neighmax, natom), molpos(3,natom), chg(natom) )
   !
   id = 0
   flag = 1
   last_CM = 0
   last_CT = 0
   first_CM = 0
   first_CT = 0
   !
   ! First H
   do ii = 1, 1
      do i = 1, 1
         id = id + 1
         element(id)   = HCHC_e(i)
         chg(id)       = HCHC_q(i)
         neigh(:, id)  = HCHC_n(:, i)
         molpos(:, id) = HCHC_p(:, i)
      enddo
      refpos = HCHC_p(:, i)
      dR = HCHC_R(1)
      nR = HCHC_R(2)
      !
   enddo
   !
   ! CH2-CH2
   do ii = 1, nx
      href  = id
      dL = CTCT_L(1) + href
      nL = CTCT_L(2)
      neigh(nR, dR) = dL
      !
      if (flag == 1) then
         first_CT = id + 1
         flag = 0
      endif
      !
      do i = 1, 6
         id = id + 1
         element(id) = CTCT_e(i)
         chg(id)     = CTCT_q(i)
         where (CTCT_n(:, i).gt.0) neigh(:, id) = href + CTCT_n(:, i)
         where (CTCT_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = refpos + CTCT_p(:, i)
      enddo
      last_CT = id - 2
      !
      neigh(nL, dL)   = dR
      dR = CTCT_R(1) + href
      nR = CTCT_R(2)
      refpos = refpos + CTCT_p(:, i)
   enddo
   !
   ! (CH=CH)m-(CH2CH2)n-
   do jj = 1, ny
      ! CH=CH
      do ii = 1, m
         href  = id
         dL = CMCM_L(1) + href
         nL = CMCM_L(2)
         neigh(nR, dR) = dL
         !
         if (flag == 1) then
            first_CM = id + 1
            flag = 0
         endif
         !
         do i = 1, 4
            id = id + 1
            element(id) = CMCM_e(i)
            chg(id)     = CMCM_q(i)
            where (CMCM_n(:, i).gt.0) neigh(:, id) = href + CMCM_n(:, i)
            where (CMCM_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = refpos + CMCM_p(:, i)
         enddo
         last_CM = id - 1
         !
         neigh(nL, dL)   = dR
         dR = CMCM_R(1) + href
         nR = CMCM_R(2)
         refpos = refpos + CMCM_p(:, i)
      enddo
      !
      ! CH2CH2
      do ii = 1, n
         href  = id
         dL = CTCT_L(1) + href
         nL = CTCT_L(2)
         neigh(nR, dR) = dL
         !
         if (flag == 1) then
            first_CT = id + 1
            flag = 0
         endif
         !
         do i = 1, 6
            id = id + 1
            element(id) = CTCT_e(i)
            chg(id)     = CTCT_q(i)
            where (CTCT_n(:, i).gt.0) neigh(:, id) = href + CTCT_n(:, i)
            where (CTCT_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = refpos + CTCT_p(:, i)
         enddo
         last_CM = 0 ! last C is not CM
         last_CT = id - 2
         !
         neigh(nL, dL)   = dR
         dR = CTCT_R(1) + href
         nR = CTCT_R(2)
         refpos = refpos + CTCT_p(:, i)
      enddo
   enddo
   !
   ! CH=CH
   do ii = 1, nz
      href  = id
      dL = CMCM_L(1) + href
      nL = CMCM_L(2)
      neigh(nR, dR) = dL
      !
      if (flag == 1) then
         first_CM = id + 1
         flag = 0
      endif
      !
      do i = 1, 4
         id = id + 1
         element(id) = CMCM_e(i)
         chg(id)     = CMCM_q(i)
         where (CMCM_n(:, i).gt.0) neigh(:, id) = href + CMCM_n(:, i)
         where (CMCM_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = refpos + CMCM_p(:, i)
      enddo
      last_CM = id - 1
      !
      neigh(nL, dL)   = dR
      dR = CMCM_R(1) + href
      nR = CMCM_R(2)
      refpos = refpos + CMCM_p(:, i)
   enddo
   !
   ! CH2
   do ii = 1, nq
      href  = id
      dL = CTHC_L(1) + href
      nL = CTHC_L(2)
      neigh(nR, dR) = dL
      !
      if (flag == 1) then
         first_CT = id + 1
         flag = 0
      endif
      !
      do i = 1, 3
         id = id + 1
         element(id) = CTHC_e(i)
         chg(id)     = CTHC_q(i)
         where (CTHC_n(:, i).gt.0) neigh(:, id) = href + CTHC_n(:, i)
         where (CTHC_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = refpos + CTHC_p(:, i)
      enddo
      last_CM = 0 ! last C is not CM
      last_CT = id - 2
      !
      neigh(nL, dL)   = dR
      dR = CTHC_R(1) + href
      nR = CTHC_R(2)
      refpos = refpos + CTHC_p(:, i)
   enddo
   !
   ! Last H
   do ii = 1, 1
      href  = id
      dL = HCHC_L(1) + href
      nL = HCHC_L(2)
      neigh(nR, dR) = dL
      !
      do i = 1, 1
         id = id + 1
         element(id) = HCHC_e(i)
         chg(id)     = HCHC_q(i)
         where (HCHC_n(:, i).gt.0) neigh(:, id) = href + HCHC_n(:, i)
         where (HCHC_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = refpos + HCHC_p(:, 3)
      enddo
      !
      neigh(nL, dL)   = dR
      dR = HCHC_R(1) + href
      nR = HCHC_R(2)
      refpos = refpos + HCHC_p(:, i)
   enddo
   ! Extra charges on first and last C 
   if (first_CT > 0) then
      chg(first_CT) = chg(first_CT) - HCHC_q(1)
   endif
   if (last_CT > 0) then
      chg(last_CT) = chg(last_CT) - HCHC_q(1)
   endif
   if (first_CM > 0) then
      chg(1) = HCHC_q(2)
      chg(first_CM) = chg(first_CM) - chg(1)
   endif
   if (last_CM > 0) then
      chg(id) = HCHC_q(2)
      chg(last_CM) = chg(last_CM) - chg(id)
   endif
   if (natom == 2) then
      chg(1) = 0.D0
      chg(2) = 0.D0
   endif
   !-------------------------------------------------------------------------
end subroutine hydrocarbon_init
